set.seed(1337)

library(VGAM)
library(gtools)
library(snow)
library(Rmpi)
library(minqa)

source("icsw.R")

###

dgp <- function(Z,X1,X2,X3,a,b,c,d,e,f,g,h) {
	N <- length(Z)
	AC <- (2 + a + b*X1 + c*X2 + d*X3 + rnorm(N,0,3) > 0)
	A <- (-2 + e + f*X1 + g*X2 + h*X3 + rnorm(N,0,3) > 0)*AC
	C <- AC - A
	
	D <- C*Z + A
	Y0 <- 5*AC + rnorm(N)
	tau<- 5*(1+X1+X2+X3)
	Y <- Y0 + tau*(D)
	
	return(list(D,Y,tau,C))
	}

iteration <- function(N) {

numbs <- 500
alpha <- 0.275

Z <- sample(c(rep(1,N/2),rep(0,N/2)))

X1 <- runif(N,-1,1)
X2 <- runif(N,-1,1)
X3 <- runif(N,-1,1)

a <- runif(1,-2,2)
b <- runif(1,-2,2)
c <- runif(1,-2,2)
d <- runif(1,-2,2)
e <- runif(1,-2,2)
f <- runif(1,-2,2)
g <- runif(1,-2,2)
h <- runif(1,-2,2)

output <- rep(NA,7)

# DGP1: Correct Model

dgpout <- dgp(Z,X1,X2,X3,a,b,c,d,e,f,g,h)

D <- dgpout[[1]]
Y <- dgpout[[2]]

Xmat <- W <- cbind(X1,X2,X3)

cscoreout <- complier.score(genoud=FALSE,D,Z,Xmat)
cscore <- cscoreout$C.score
cscore[cscore < .Machine$double.eps^.5] <- .Machine$double.eps^.5
qcscore <- quantile(cscore,1/N^alpha)
cscore[cscore < qcscore] <- qcscore

output[1] <- mean(dgpout[[3]])
output[2] <- lm(Y~Z+Xmat)$coefficients["Z"]
output[3] <- lm(Y~D+Xmat)$coefficients["D"]
output[4] <- tsls.wfit(cbind(1,D,Xmat),Y,cbind(1,Z,Xmat),weights=rep(1,N))$coefficients["D"]
output[5] <- tsls.wfit(cbind(1,D,Xmat),Y,cbind(1,Z,Xmat),weights=1/cscore)$coefficients["D"]
output[6] <- mean(dgpout[[4]])
output[7] <- N

BScoefs <- rep(-999,numbs) # change back

for (i in 1:numbs) {
	bssamp <- sample(1:N,replace=TRUE)
	cscoreBS <- complier.score(genoud=FALSE,D[bssamp],Z[bssamp],Xmat[bssamp,])$C.score
	cscoreBS[cscoreBS < .Machine$double.eps^.5] <- .Machine$double.eps^.5
	qcscoreBS <- quantile(cscoreBS,1/N^alpha)
cscoreBS[cscoreBS < qcscoreBS] <- qcscoreBS
	BScoefs[i] <- tsls.wfit(cbind(1,D,Xmat)[bssamp,],Y[bssamp],cbind(1,Z,Xmat)[bssamp,],weights=1/cscoreBS)$coefficients["D"]
	}

c(output,quantile(BScoefs,c(0.025,0.05,0.95,0.975)))

return(c(output,quantile(BScoefs,c(0.025,0.05,0.95,0.975))))
}

cl <- makeSOCKcluster(7)

clusterEvalQ(cl,source("icsw.R"))
clusterEvalQ(cl,library(minqa))
clusterEvalQ(cl,library(VGAM))
clusterEvalQ(cl,library(gtools))

clusterExport(cl,ls())

numiter <- 500

Nlist <- c(500, 2500, 5000)

clusterSetupRNG(cl)

simoutput <- clusterApplyLB(cl,sort(rep(Nlist,numiter)),iteration)

save.image("Bootstrap-Sims.rData")

simoutput <- do.call(cbind,simoutput)

Nno <- simoutput[7,] == 500

# 90% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[9,Nno] & simoutput[1,Nno] <= simoutput[10,Nno])

# 95% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[8,Nno] & simoutput[1,Nno] <= simoutput[11,Nno])

Nno <- simoutput[7,] == 2500

# 90% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[9,Nno] & simoutput[1,Nno] <= simoutput[10,Nno])

# 95% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[8,Nno] & simoutput[1,Nno] <= simoutput[11,Nno])

Nno <- simoutput[7,] == 5000

# 90% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[9,Nno] & simoutput[1,Nno] <= simoutput[10,Nno])

# 95% Interval Coverage
mean(simoutput[1,Nno] >= simoutput[8,Nno] & simoutput[1,Nno] <= simoutput[11,Nno])

stopCluster(cl)
